#include "../Iteration_Params.hxx"

Okada::Displacement Iteration_Params::uc
(const double &z,
 const FTensor::Tensor1<double,3> &dislocation,
 const double &alp4, const double &alp5,
 const double &sin_dip, const double &cos_dip,
 const double &cos_dip2,
 const double &sin_dip2,
 const double &sin_cos_dip) const
{
  const double c(d+z), x53((8*r2+9*r*xi+3*xi2)*x11*x11*x11/r2),
    y53((8*r2+9*r*et+3*et2)*y11*y11*y11/r2),
    h(q*cos_dip-z), z32(sin_dip/r3-h*y32), z53(3*sin_dip/r5-h*y53), y0(y11-xi2*y32),
    z0(z32-xi2*z53), ppy(cos_dip/r3+q*y32*sin_dip), ppz(sin_dip/r3-q*y32*cos_dip),
    qq(z*y32+z32+z0), qqy(3*c*d/r5-qq*sin_dip), qqz(3*c*y/r5-qq*cos_dip+q*y32),
    xy(xi*y11), qy(q*y11), qr(3*q/r5),
    cdr((c+d)/r3), yy0(y/r3-y0*cos_dip);

  Okada::Displacement result;
  FTensor::Index<'a',3> a;
  FTensor::Index<'b',3> b;
  result.first(a)=0;
  result.second(a,b)=0;

  const double pi(4*atan(1.0));
  FTensor::Tensor1<double,3> first;
  FTensor::Tensor2<double,3,3> second;
  /* STRIKE-SLIP CONTRIBUTION */
  if(dislocation(0)!=0)
    {
      first(0)= alp4*xy*cos_dip           -alp5*xi*q*z32;
      first(1)= alp4*(cos_dip/r+2*qy*sin_dip) -alp5*c*q/r3;
      first(2)= alp4*qy*cos_dip           -alp5*(c*et/r3-z*y11+xi2*z32);
      second(0,0)= alp4*y0*cos_dip                  -alp5*q*z0;
      second(0,1)=-alp4*xi*(cos_dip/r3+2*q*y32*sin_dip) +alp5*c*xi*qr;
      second(0,2)=-alp4*xi*q*y32*cos_dip            +alp5*xi*(3*c*et/r5-qq);
      second(1,0)=-alp4*xi*ppy*cos_dip    -alp5*xi*qqy;
      second(1,1)= alp4*2*(d/r3-y0*sin_dip)*sin_dip-y/r3*cos_dip
        -alp5*(cdr*sin_dip-et/r3-c*y*qr);
      second(1,2)=-alp4*q/r3+yy0*sin_dip  +alp5*(cdr*cos_dip+c*d*qr-(y0*cos_dip+q*z0)*sin_dip);
      second(2,0)= alp4*xi*ppz*cos_dip    -alp5*xi*qqz;
      second(2,1)= alp4*2*(y/r3-y0*cos_dip)*sin_dip+d/r3*cos_dip -alp5*(cdr*cos_dip+c*d*qr);
      second(2,2)=         yy0*cos_dip    -alp5*(cdr*sin_dip-c*y*qr-y0*sin_dip2+q*z0*cos_dip);

      result.first(a)=first(a)*dislocation(0)/(2*pi);
      result.second(a,b)=second(a,b)*dislocation(0)/(2*pi);
    }

  /* DIP-SLIP CONTRIBUTION */
  if(dislocation(1)!=0)
    {
      first(0)= alp4*cos_dip/r -qy*sin_dip -alp5*c*q/r3;
      first(1)= alp4*y*x11       -alp5*c*et*q*x32;
      first(2)=     -d*x11-xy*sin_dip -alp5*c*(x11-q2*x32);
      second(0,0)=-alp4*xi/r3*cos_dip +alp5*c*xi*qr +xi*q*y32*sin_dip;
      second(0,1)=-alp4*y/r3     +alp5*c*et*qr;
      second(0,2)=    d/r3-y0*sin_dip +alp5*c/r3*(1-3*q2/r2);
      second(1,0)=-alp4*et/r3+y0*sin_dip2 -alp5*(cdr*sin_dip-c*y*qr);
      second(1,1)= alp4*(x11-y*y*x32) -alp5*c*((d+2*q*cos_dip)*x32-y*et*q*x53);
      second(1,2)=  xi*ppy*sin_dip+y*d*x32 +alp5*c*((y+2*q*sin_dip)*x32-y*q2*x53);
      second(2,0)=      -q/r3+y0*sin_cos_dip -alp5*(cdr*cos_dip+c*d*qr);
      second(2,1)= alp4*y*d*x32       -alp5*c*((y-2*q*sin_dip)*x32+d*et*q*x53);
      second(2,2)=-xi*ppz*sin_dip+x11-d*d*x32-alp5*c*((d-2*q*cos_dip)*x32-d*q2*x53);

      result.first(a)+=first(a)*dislocation(1)/(2*pi);
      result.second(a,b)+=second(a,b)*dislocation(1)/(2*pi);
    }

  /* TENSILE-FAULT CONTRIBUTION */
  if(dislocation(2)!=0)
    {
      first(0)=-alp4*(sin_dip/r+qy*cos_dip)   -alp5*(z*y11-q2*z32);
      first(1)= alp4*2*xy*sin_dip+d*x11 -alp5*c*(x11-q2*x32);
      first(2)= alp4*(y*x11+xy*cos_dip)  +alp5*q*(c*et*x32+xi*z32);
      second(0,0)= alp4*xi/r3*sin_dip+xi*q*y32*cos_dip+alp5*xi*(3*c*et/r5-2*z32-z0);
      second(0,1)= alp4*2*y0*sin_dip-d/r3 +alp5*c/r3*(1-3*q2/r2);
      second(0,2)=-alp4*yy0           -alp5*(c*et*qr-q*z0);
      second(1,0)= alp4*(q/r3+y0*sin_cos_dip)   +alp5*(z/r3*cos_dip+c*d*qr-q*z0*sin_dip);
      second(1,1)=-alp4*2*xi*ppy*sin_dip-y*d*x32
        +alp5*c*((y+2*q*sin_dip)*x32-y*q2*x53);
      second(1,2)=-alp4*(xi*ppy*cos_dip-x11+y*y*x32)
        +alp5*(c*((d+2*q*cos_dip)*x32-y*et*q*x53)+xi*qqy);
      second(2,0)=  -et/r3+y0*cos_dip2 -alp5*(z/r3*sin_dip-c*y*qr-y0*sin_dip2+q*z0*cos_dip);
      second(2,1)= alp4*2*xi*ppz*sin_dip-x11+d*d*x32
        -alp5*c*((d-2*q*cos_dip)*x32-d*q2*x53);
      second(2,2)= alp4*(xi*ppz*cos_dip+y*d*x32)
        +alp5*(c*((y-2*q*sin_dip)*x32+d*et*q*x53)+xi*qqz);

      result.first(a)+=first(a)*dislocation(2)/(2*pi);
      result.second(a,b)+=second(a,b)*dislocation(2)/(2*pi);
    }
  return result;
}
